SUBROUTINE ZERODIM_MODEL &
& (h, dt, &
& shortwave, longwave, tempair, humair, pressure, &
& uwind, vwind, zref, Temp0)

! Subroutine T_0DIM calculates the lake temperature
! on the new timestep according to zero-dimensional model
! (well mixed lake)

use PHYS_CONSTANTS, only : &
& cw_m_row0, &
& albedoofwater, &
& emissivityofwater, &
& sigma, &
& roughness0, &
& Kelvin0, &
& pref, &
& Rd, Rd_d_cp, Rd_d_Rwv, &
& esatsurf0, &
& aMagw, bMagw, &
& Lwv, cp

implicit none

! Input/output variables

real(8), intent(in) :: h ! Lake depth, m
real(8), intent(in) :: dt ! Timestep, s
real(8), intent(in) :: shortwave ! Shortwave radiation flux, W/m**2
real(8), intent(in) :: longwave ! Longwave radiation flux, W/m**2
real(8), intent(in) :: tempair ! Air temperature, Celsius
real(8), intent(in) :: humair ! Specific humidity, kg/kg
real(8), intent(in) :: pressure ! Air pressure, Pa
real(8), intent(in) :: uwind ! X-wind component, m/s
real(8), intent(in) :: vwind ! Y-wind component, m/s
real(8), intent(in) :: zref ! The height of measurements 
                            ! (or of the lowermost atmospheric model level), m
real(8), intent(inout) :: Temp0 ! Lake temperature, Celsius

! Local variables and parameters
real(8), parameter :: minwind = 1.d-2
real(8), parameter :: relhums = 1.
integer(4), parameter :: itdrag = 10

real(8) :: eflux, hflux, Elatent, surfrad
real(8) :: wr, TET2, TET1, ro, esatsurf, humsurf
real(8) :: c_u, c_t
real(8) :: bx(7), bix(11)

wr       = dmax1(dsqrt(uwind**2 + vwind**2), minwind)
TET2     = (tempair + Kelvin0)*(pref/pressure)**Rd_d_cp
TET1     = (Temp0 + Kelvin0)*(pref/pressure)**Rd_d_cp  
ro       = pressure/(Rd*(tempair + Kelvin0))  
esatsurf = esatsurf0*10.**(aMagw*Temp0/(bMagw + Temp0))
humsurf  = Rd_d_Rwv/pressure*esatsurf*relhums

! Businger-Dayer, Beljaars parameterization for exchange coefficients
bx(1) = wr
bx(2) = TET2
bx(3) = TET1
bx(4) = humair 
bx(5) = humsurf
bx(6) = zref
bx(7) = roughness0
call DRAGVL (bx, bix, itdrag)
c_u = bix(8)
c_t = bix(9)
hflux   = - cp*ro*c_u*c_t*wr*(TET2 - TET1)
Elatent = - Lwv*ro*c_u*c_t*wr*(humair - humsurf) 
surfrad = sigma*emissivityofwater*(Temp0 + Kelvin0)**4

eflux = shortwave*(1. - albedoofwater) + longwave - surfrad - hflux - Elatent
Temp0 = Temp0 + dt*eflux/(cw_m_row0*h)

END SUBROUTINE ZERODIM_MODEL